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ABSTRACT 

Based  on  the  original  development  of  radial  mixing  effects  by  Gallimore,  the  development 
of  the  governing  equations  for  hub-to-tip  streamline  curvature  solutions  in  tubomachinery 
passages  and  blade  rows  and  incorporating  viscous  radial  mixing  is  presented.  Arbitrary 
passage  geometries  and  quasi-orthogonal  descriptions  are  accounted  for.  FORTRAN 
program  coding  is  carried  out  for  the  special  case  of  unbladed,  axial-flow  passages.  Flow 
solutions  are  obtained  by  iteration  based  on  fixed  streamline  patterns  in  the  meridional 
plane  followed  by  further  iterations  involving  readjustment  of  the  streamlines  until  overall 
convergence  is  obtained.  Simple  examples  are  presented  to  illustrate  the  solution 
technique.  This  work  was  carried  out  under  AFOSR  Grant  F49620-92-J-(X)26. 
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SYMBOLS  AND  NOTATION 

A,  B  Parameter  used  to  integrate  the  momentum  equation 

B  unit  vector  tangent  to  the  mean  blade  surface  and  which  lies  in  the  q-0 

plane  at  an  angle  e  from  the  q-o  direction 
E„,  Ee,  E„  viscous  body  force  components 
f.,  f®,  fm  distributed  blade  body  force  components 
fp  force  vector  normal  to  the  blade  surface 

g  gravitational  acceleration 

H  stagnation  enthalpy 

h  static  enthalpy 

i  quasi-orthogonal  index 

j  streamline  index 

k  heat  transfer  coefficient 

m  meridional  direction 

m  mass  flow  rate 

n  normal  direction 

p  static  pressure 

r  radial  direction 

T°  stagnation  pressure 

T  static  temperature 

z  axial  direction 

Vi,  V„,  Vg,  V,  velocities 


Sc  Schmidt  number  Sc  =  ^/(pc) 

Pr  Prandtl  number  Pr  =  p  Cp/k 

e„  6,,  unit  vectors 

e„,  Cgg,  e™  normal  strain  rates 

Cna,  6,16,  Girt  shear  strain  rates 

r^  radius  of  curvature 

q  heat  flux 

s  entropy 

GREEK 

B  angle  measured  from  the  meridional  direction  to  a  unit  vector  tangent  to 

the  turbomachine  blade  (see  Fig.  Al) 

Y  quasi-orthogonal  lean  angle  measured  from  the  radial  direction 


a.  Q. 


V 


e  mixing  coefficient 

0  tangential  direction 

viscosity 
density 

(^ee,  <^mm  stess  tensor  comonenets 

angle  measured  from  the  axial  direction  to  a  line  tangent  to  a  streamline 
viscous  dissipation  function 
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1.  INTRODUCTION 

For  many  ye wrs,  axisymmetric  solutions  obtained  from  streamline  curvature 
analyses  have  been  a  standard  approach  in  analyzing  axial  compressor  or  turbine  design 
and  performance.  In  this  approach,  the  flow  calculations  are  reduced  as  a  practical 
consideration  to  a  mathematically  two-dimensional  problem  in  which  the  through-flow  is 
assumed  to  be  steady,  inviscid,  and  confined  to  concentric  stream  surfaces.  As  a  result, 
no  transport  of  mass,  momentum,  or  energy  takes  place  across  these  concentric  surfaces, 
and  thus  no  accounting  of  essentially  three-dimensional  effects  is  made  due  to  secondary 
flows,  blade  or  endwall  boundary  layer  flows,  or  tip  leakage  flows.  In  particular, 
comparisons  with  detailed  interstage  traverses  have  shown  that  calculated  results  do  not 
accurately  model  radial  profiles  of  stagnation  temperature,  pressure,  and  flow  angle  as  the 
profiles  develop  through  the  machine.  Thus  inclusion  of  spanwise  mixing  terms  in  the 
through-flow  calculations  would  serve  as  important  function  in  improving  the  modeling  of 
radial  profiles  and  allow  more  realistic  modeling  of  losses.  Streamline  curvature 
calculations  are  carried  out  in  an  iterative  fashion  to  satisfy  momentum  and  continuity 
requirements  along  quasi-oithogonals  established  in  the  two-dimensional  field.  These 
calculations  have  proved  to  be  acceptable  in  terms  of  stability  and  convergence,  at  least  in 
the  case  of  hub-to-tip  solutions  conducted  in  the  meridional  plane.  This  quasi  3-D 
approach  is  very  likely  to  continue  as  the  standard  method  of  aerodynamic  analysis  in 
turbomachinery  for  some  years  to  come,  although  fully  3-D  codes  are  used  on  occasion  to 
further  refine  details  of  the  flow. 
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Practical  approaches  to  account  for  secondary  flow  effects  have  had  a  considerably 
heuristic  goal  in  devising  reasonably  simple  but  adequate  models.  Wennerstrom  (1991),  in 
reviewing  modeling  ideas  that  have  been  proposed  over  the  last  20  or  25  years,  points  out 
that  an  approach  based  on  transport  phenomena  has  been  used  to  superimpose  semi- 
empirical  models  of  secondary  flow  features  as  corrections  to  otherwise  classical 
axisymmetric  solutions.  Wennerstrom  has  termed  these  modeling  attempts  as  hybrid 
solutions  since  they  are  not  truly  two-dimensional,  nor  are  they  fully  three-dimensional. 

By  necessity,  the  objective  is  not  to  deal  with  the  three-dimensional  details  of  the  flow, 
but  to  deal  with  the  overall  computation  for  arbitrary  geometries.  The  major  contributors 
to  the  modeling  efforts  in  recent  years  have  been  Adkins  and  Smith  (1982),  who  proposed 
that  secondary  flows  were  responsible  for  spanwise  mixing,  and  Gallimore  and  Cumpsty 
(1986),  who  viewed  turbulent  diffusion,  rather  than  secondary  flow,  to  be  the  dominant 
mechanism  of  spanwise  mixing.  Leylek  and  Wisler  (1991)  have,  through  extensive 
numerical  analysis  and  experiments,  largely  resolved  the  issue  of  the  dominant  mixing 
mechanism.  They  concluded  that  the  relative  importance  of  secondary  flow  and  turbulent 
diffusion  depend  -non  stage  loading,  solidity,  and  aspect  ratio.  Thus,  it  would  appear 
that  both  mechanisms  should  be  incorporated  in  the  performance  calculation  process.  It  is 
worth  noting  that  if  one  is  able  to  predict  with  reasonable  assurance  the  performance  of  a 
given  machine,  then  one  has  the  most  powerful  design  tools  available.  Howard  and 
Gallimore  (1992)  extended  the  turbulent  diffusion  model  of  Gallimore  and  Cumpsty 
(1986)  to  include  a  viscous  wall. 
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The  objective  of  the  research  described  in  this  work  was  to  develop  the  governing 
equations  for  hub-to-tip  flow  analysis  in  intrinsic  coordinates  using  the  streamline 
curvature  throughflow  method  incorporating  the  turbulent  diffusion  model  of  spanwise 
mixing  developed  by  Gallimore  and  Cumpsty. 

A  special  case  of  the  governing  equations  for  unbladed,  axial-flow  passages  was 
then  used  to  create  a  FORTRAN  computer  program  to  demonstrate  the  iterative  solution 
technique. 
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2.  GOVERNING  EQUATIONS  FOR  HUB-TO-TIP  FLOW  ANALYSIS 

2.1  The  Momentum  Equations 

Turbomachinery  flow  is  viscous,  compressible,  unsteady,  and  highly  three- 
dimensional.  One  method  often  used  in  turbomachinery  design  and  analysis  is  to  assume 
steady,  axisymmetric  flow  through  the  annular  passage,  allowing  the  flow  to  be  treated  in 
terms  of  a  so-called  hub-to-tip  analysis  in  the  meridional  plane.  Thus  the  flow  analysis  is 
reduced  to  one  involving  only  two  space  dimensions.  In  the  case  of  bladed  regions,  these 
flow  simplifications  are  equivalent  to  assuming  an  infinite  number  of  infinitely  thin 
blades. 

Figure  2.1  shows  the  meridional  plane;  a  cylindrical  coordinate  system  (r,  0,  z)  is 
indicated  with  the  z  axis  aligned  along  the  axis  of  the  machine,  and  r  pointing  in  the 
radial  direction.  The  0  or  tangential  direction  points  out  of  the  paper  normal  to  the 
meridional  plane.  Also  shown  in  Figure  2.1  is  an  intrinsic  coordinate  system  (n,  0,  m) 
with  the  m  and  n  coordinates  lying  in  the  plane.  The  m  coordinate  corresponds  to  the 
meridional  streamlines  in  the  flow,  and  the  n  coordinate  is  orthogonal  to  them.  Thus  the 
flow  velocity  can  be  expressed  as 


(2.1) 


where  V,  is  the  meridional  velocity,  and  V,  is  the  tangential  velocity.  The  angle  <f)  is  the 
slope  of  the  streamlines  measured  from  the  z-axis,  and  e,  and  are  unit  vectors  along  the 
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m  and  n  coordinate  directions.  The  («,  0,  m)  coordinates  will  be  used  in  the  development 
of  the  analysis  which  follows. 

In  Figure  2.2,  the  meridional  plane  is  shown  again  with  representative  hub  and  tip 
boundaries  of  the  annulus,  and  with  rotor  and  stator  blade  rows  included.  In  intrinsic 
coordinates,  the  normal,  tangential,  and  meridional  momentum  equations  may  be  written 
(Tsien,  1958)  as 


cos<|)  -  V. 


2^ 

* 

cm 


do 

nn 

+  0 

'sin4>  ^ 

d^" 

dm 

dn 
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1  r 
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P/e 


(2.3) 


7 


(2.4) 


where  p  is  density,  the  various  o  terms  are  the  stress  tensor  components,  the  /terms  are 
distributed  blade  body  forces,  and  V®,  and  r  are  as  previously  defined.  See 
Appendix  A  for  the  derivation  of  the  distributed  blade  body  force  components. 

The  stress-strain  rate  relationships  for  a  compressible  fluid  may  be  found  in  White 
(1991).  Using  thermodynamic  pressure  and  viscosity,  the  stress-strain  rate  relationships 
written  in  terms  of  the  intrinsic  coordinates  (n,  0,  m)  are: 


=  ~p 

®0e 

=  -p 

^mm 

=  ~p 

+  2iie^ 

=  2pe^ 

=  Oe, 

=  2iie^ 

®«e 

II 

Q 

? 

=  2ne^ 

The  strain  rates  may  be  expressed  as  functions  of  the  velocity  gradients  (White,  1991). 
Substituting  for  the  strain  rates  in  Eqn.  (2.5),  the  stress-strain  rate  equations  take  the 


form: 
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°nn=  -P  ^ 
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“  •‘U 


dn  dm) 

—  C0S<|) 
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(2.6) 


o«0  ~  Oe«  ~ 


dV,  V,  . 

- -  — sm<b 

dm  r 


These  expressions  for  the  stresses  may  be  substituted  into  Eqns.  (2.2)  -  (2.4),  leaving  the 
equations  of  motion  as  functions  of  p,  p,  4>,  fi,  V„  and  V#,  and and  f„. 

As  an  approach  to  viscous  analysis  of  hub-to-tip  flows,  Gallimore  (1985)  proposed 
a  spanwise  mixing  model  in  cylindrical  coordinates  of  the  velocity  and  temperature  fields 
by  retaining  in  the  momentum  equations  the  normal  stresses,  a„,  a^,  Oa,  and  only  the 
shear  stresses  acting  on  the  inner  and  outer  cylindrical  surfaces  of  the  fluid  element, 
and  a„.  He  argued  also  that,  because  mixing  is  assumed  to  occur  across  the  cylindrical 
surfaces,  the  normal  strain  rates  in  the  axial  and  tangential  directions  (e^  and  Cee)  could 
be  neglected.  Heat  transfer  due  to  mixing  would  be  modeled  by  retaining  only  radial 
temperature  gradients  in  the  energy  equation. 

A  spanwise  mixing  model  for  intrinsic  coordinates  may  be  derived  in  a  similar 
fashion  with  mixing  assumed  to  occur  across  the  meridional  stream  surfaces.  Therefore, 
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the  shear  stresses  ct.#  and  ff„„  are  retained  while  the  remaining  shear  stresses  are  neglected 
and  the  normal  stresses,  <t„„,  and  a,*,  are  retained  with  the  strain  rates  e„„  and  eee 
neglected.  Heat  transfer  due  to  mixing  is  modeled  by  retaining  only  normal  temperature 
gradients  in  the  energy  eQuation.  Figure  2.3  shows  a  fluid  element  in  intrinsic  coordinates 
with  the  stress  terms  retained  by  the  spanwise  mixing  model.  The  momentum  equations 
obtained  from  Eqns.  (2.2)  -  (2.4)  and  (2.6)  incorporating  the  spanwise  mixing  model  in 
intrinsic  coordinates  are; 


[^e  . 

=  t-  2- 

dn) 

-  V  .^1 

- C08<P  - 

r 

-  Y  — 

— 

om 

dn 

dn 

,  dn 
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(2.9) 


The  notation  in  these  equations  may  be  simplified  to 


"•dm  r  ^  ^ 


p  dn 


(2.10) 


y.  d{rV,) 

r  dm 


~  ^9  "^/e 


(2.11) 


dV„  F  2  ,  g_ 

ym-^  -  — sm<|>  =  -1^  +  +  A 


dm  r 


(2.12) 


where  the  E„  £#,  and  E„  terms  represent  viscous  body  forces  per  unit  mass  and  are 


written  as: 
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COS<|)  _  0<t) 
r  dm 


(2.13) 


2.2  The  Radial  Equilibrium  Equation 

The  determination  of  flow  variables  in  the  meridional  plane  is  accomplished  by  the 
iterative  solution  of  the  uncoupled  momentum  and  continuity  equations.  The  solution  of 
the  normal  and  meridional  momentum  equations  is  complicated  by  the  fact  that  both  the 
n  and  m  directions  change  throughout  the  flow  field  as  the  computation  approaches 
convergence.  The  flow  field  solution  may  be  simplified  by  constructing  fixed  quasi- 
orthogonal  {q-o)  coordinates  approximately  perpendicular  to  the  streamlines.  The  n 
direction  may  then  be  eliminated  by  transforming  the  momentum  equations  to  the  m  and 
q-o  directions.  A  flow  field  grid  is  then  formed  by  the  streamlines  and  the  quasi- 
orthogonals.  Figure  2.4  shows  the  coordinate  systems  in  the  meridional  plane,  including 
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the  defined  q-o.  The  angle,  y,  between  the  r  and  q-o  directions  is  taken  to  be  positive 
when  measured  as  shown. 

The  chain  rule  for  a  function  of  two  variables  may  be  written  to  express  the 
directional  derivative  along  the  q-o  as: 

d( )  ^  ^  ^2  14) 

dq  dn  dq  dm  dq 

In  addition,  the  following  geometric  relationships  can  be  found  from  Figure  2.4: 

=  co8(4)  +  y)  =  sin(<|)  +  y)  (2.15) 

dq  dq 

Combining  Eqns.  (2.14)  and  (2.15)  and  solving  for  d(  )/dn  yields 

^  =  sec(<|>  +  Y)-^  -  (2-16) 

dn  dq  dm 

This  relationship  is  used  in  the  development  below  to  transform  the  momentum 
component  equations  to  the  (m,  q)  coordinate  system. 

From  Gibbs’  Tds  relationship  we  can  write 

^  =  -  Tds  +  dff  -  VdV  (2-17) 

P 


where  H  is  the  stagnation  enthalpy  and  V  is  the  flow  velocity.  The  pressure  gradients  in 
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the  m  and  n  directions  therefore  are 


_1.JP 

_  _jds 

.  ^  -  y^y 

p  dm 

dm 

dm  dm 

1  dp 

=  -T^ 

p  dn 

dn 

dn  dn 

(2.18) 


Substituting  Eqns  (2.18)  into  the  normal  and  meridional  momentum  equations,  Eqns. 
(2.10)  and  (2.12),  the  momentum  equations  become: 
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Using  Eqn.  (2.16)  to  replace  the  n  directional  gradients  in  the  normal  momentum 
equation  with  gradients  in  the  m  and  q  directions,  we  obtain 


y  2  ^ 
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.dH  ,  .dH 

=  -sec(<j>  +  Y)—  +  tan(<|)  +  Y)— 
dq  dm 
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The  m  and  n  momentum  equations  may  now  be  combined  by  multiplying  the  m 
momentum  equation  by  tan(<^  +  y)  and  adding  the  result  to  the  n  momentum  equation. 
Further  rearrangement  gives 


tan(<|)  +  y) 


V  — -  -  — siiKb  -  E  -  f 

m  dm  r  ” 


+  V^—  -  — cos4)  = 
"aw  r 

dV. 


-sec(<|>  +  Y)-^  +  sec(^+y)T^  +  sec(<|>  +  Y)t^e-^ 


as 


(2.21) 


dq 


dq 


dq 


dV 

sec(<|)+Y)t^„-^ 


Multiplying  Eqn.  (2.21)  by  cos(</>  +  y)  and  rearranging  gives: 


dq  dm  dm  dq  dq 


dV^ 


sin(<t)  +  y){£„  +/«)  -  cos(<|)  +  y){E„  +/„) 


(2.22) 


^9 


dVa 


r  vr 

-  cos(<|>  +  y)cos<|)—  -  sin(4)  +  Y)sm<l)—  -  K«-^ 
r  r  dq 


The  radius  of  curvature,  r„  of  the  meridional  streamline  is  defined  by: 


1  = 

dm 


(2.23) 


Using  this  result  in  Eqn  (2.22)  along  with  further  rearrangement,  and  recognizing  that 
dr/dq  =  cos  y  >  we  get  the  final  form  of  radial  equilibrium  equation 
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dV  dV  V  ^  Au 

V„-^  ^  cos«|)n)^  ^  ^  -  Tf 

dq  dm  dq  dq 

-  sm^y)[E„^f„]  -  cos(<t)+Y)p„+/„]  - 


(2.24) 


Equation  (2.16)  may  be  used  to  eliminate  the  n  directional  gradient  terms  in  the 
normal,  tangential,  and  meridional  viscous  body  force  equations,  leaving 


2sec(<l)+Y)— 

dq' 


\i  I^„|sec(<l)  +  Y)^  -  tan(<j)  +  y)-^  j 

nK^|sec((l)  +  Y)-^ -tan(<l>  +  Y)-^j  + 
sec(<|)  +  Y)-^  -tan(<|)  +  Y)-^  fsec(<|)  +  Y)^  -  tan(<j)  +  y)-^ 

C/ftt  wfft  ^  y  09t% 


2tan(<|>  +  Y)— 
dm 


(2.25) 
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^0  =  - 

sec(<|)+Y)-|- 

pr  sec(4)+Y)^ 

(v) 

’'6 

3  /  'j'] 

-tan(<|)+Y)^  — 

p 

dq 

1  dq 

i  r  J 

dm\  r  ))\ 

-  t2n(<l>+Y)^ 

owl 


V^r 


'  /A 

sec(<l)  +  Y)^  — 

r 


-tani^+y)  — 
dm 


V  r  , 


+  nr[^^--^l(sec(<|)+Y)-^[— 

\  r  om)\ 


(2.26) 


F  -  1 
P 


a  oK  dK 

sec{<t>  +  Y)—  sec(<|)  +  Y)—  -tan(<|>+Y)—  "  „ 

dq[  \  dq  dm  dm)\ 


-  tani^^y)— 

dm 


sec(<j)+Y)-^  -tan(<|>+Y)^  -  y„^ 


dq 


dm  dm)} 


+  -  ■^][sec(<|>  +  Y)-^  -  taii(<t)  +  y)^ 

\  r  dm )\  dq  dm  dm) 

_ 

-  2nK„|sec(<|>+Y)-^ 


(2.27) 


Because  momentum  and  temperature  mixing  are  assumed  to  occur  only  in  the 
normal  direction,  the  terms  in  the  viscous  body  force  equations  which  contain  velocity 
gradients  in  the  meridional  direction  will  be  neglected,  leaving 
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=  -|2sec(4)+Y)-^| 

Pi  dq 


P  K«fsec(<l>  +  Y)-^  -  tan((t)  +  y)|^ 
\  dq  dm 


-  2tan(<|)  +  Y)— 

dm 


P^«^sec(4)  +  Y)-^ -tan(<|)  +  Y)-^j  + 

n[sec(4)  +  Y)-t-^  Y)-^  -  tan(<|)  +  y)-^ 

\  dq  dm)\  dq  dm) 

+2 Ilk'  ^ISS^  - 

"awi  r  dm}] 


E  =1^ 

®  p 


sec(4)+Y)— 

dq 


-  tan(<|>+Y)^ 

dm 


P'" 


Pn 


sec(<j)  +  Y)- 


/cos<t)  _  a<t) 
\  r  dm 


^_a 

Vl 

dq 

r 

{  r 

1 

(v) 

J 

r  ) 

T 

K 

r 


|tan(<|)+Y)^ 

2  dm 


n 


M 


dr 

^tan(<l>  +  Y)-^ 

dm 


-|tan(<J)  +  Y)^ 

dm 


n 


K  =  - 

P 


sec(4>  +  Y) 


dq 


-  taii(<|>+Y)^ 
dm 


2\iV„  [sec{<J)  +  Y)-^  -  tan(4>  +  y)|^] 

V  dq  dm) 


(2.28) 


(2.29) 


(2.30) 
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In  the  iterative  solution  of  the  radial  equilibrium  equation,  the  first  two  terms  on 
the  right  hand  side  of  Eqn  (2.24)  use  values  of  V;,,  (ji  +  y,  and  r,  from  a  previous  solution 
ot  at  flow  field.  In  order  to  close  the  system,  equations  are  needed  for  the  entropy  and 
siagnition  enthalpy  gradients.  These  equations  are  derived  in  the  following  sections. 

2.3  The  Energy  Equation 

Expressed  in  terms  of  enthalpy,  the  energy  equation  may  be  written: 

-  v?  +  *  (2.31) 

Dt  Dt 

where  q  is  heat  flux  and  $  is  the  dissipation  function.  Gibbs’  Tds  relationship,  Eqn. 
(2.17),  may  be  used  to  eliminate  the  pressure  derivative  in  the  energy  equation 
yielding: 

pr—  =  -V<7  +  ®  (2.32) 

Dt 

The  left  hand  side  of  this  equation  may  be  expanded  by  using  the  definition  of  the 
substantial  derivative,  which,  due  to  the  assumption  of  steady  flow,  yields 

-  -V  ?*  «  (2.33) 

The  intrinsic  coordinate  forms  of  the  del  operator  and  the  divergence  are  needed  in 
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order  to  expand  the  heat  transfer  term.  The  del  operator  is  defined  as 


V  = 


.  a 

- 

"am 


+  e. 


_a 

an 


(2 


TI-«> 


Retaining  only  the  n  directional  temperature  gradients,  V  T  gives 


VT  = 


(2.35) 


Assuming  that  conduction  is  modeled  by  Fourier’s  Law  ,  the  gradient  of  the  heat  flux 
may  expanded  to  give 

V-(-A:)Vr  =  V{-k)'Vr  +  (-it)V-Vr  (2.36) 

where  k  is  the  coefficient  of  thermal  conductivity.  The  coordinate  system  geometry  gives 
the  following  relationships  for  the  derivatives  of  the  unit  vectors 

^  =  -^e  —  =  -^e  (2.37) 

dm  am  "  dn  an  " 


Combining  Eqn.  (2.34)  -  (2.37)  gives 


-V^  =  -V-(-itW) 


dm  dn  an\  an) 


(2.38) 


The  second  term  on  the  right  hand  side  of  the  energy  equation,  Eqn.  (2.31),  is  the 
dissipation  function,  which,  written  in  terms  of  strain  rates  (Hughs  and  Gaylord,  1964) 
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for  the  intrinsic  coordinates,  {n,  0,  m)  is 

The  dissipation  function  may  be  simplified  by  using  the  mixing  model  assumption  in 
which  the  normal  strain  rates  and  dV/dm  may  be  neglected.  With  these  assumptions,  the 
dissipation  function  becomes; 


(2.40) 


Substituting  Eqns.  (2.38)  and  (2.40),  the  energy  equation,  Eqn.  (2.31)  may  be 
written: 

*  .  ^  (2.41) 

dm  pTV^dmdn  pTV^dn\  dn)  pTV^  dm 

where  dSeldm  is  an  empirically  based  entropy  gradient  which  accounts  for  losses 
generated  by  flow  over  the  blades. 

Equation  (2.16)  is  used  to  transform  Eqn.  (2.38)  from  the  (n,  m)  to  the  {m,  q) 


coordinate  system  yielding 
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ds  _  _ 

dm  p  TV 


tan(<t> 

9  TV 


^^dql  pTV^  dq\^^  dql 

^Y)  d  { 


A:sec(<|)+Y)^| 


PTV„[[ 


sec(<|)+Y) — -  — ^sm4> 

2 

+ 

sec(4>*Y)^  - 

IL  ^  ^  J 

dq  "  dm\ 

(2.42) 


— sin<t) 
r 


dm 


2.4  The  Stagnation  Enthalpy  Equation 
The  definition  of  stagnation  enthalpy  may  be  combined  with  Gibbs’  Tds 
relationship  to  give 


j.  ^  _  dH  y  dV^  _  _  jds 

p  dm  dm  “"dm  ^  dm  dm 


(2.43) 


Substituting  Eqn.  (2.43)  into  the  meridional  momentum  equation,  Eqn.  (2.12),  yields 


oin/k  -  r*' 

a;:  ■ 

cvn  ufft 


-  Va 


dm 


+  E  +/ 

m  •fm 


(2.44) 


Solving  Eqn.  (2.44)  for  dHIdm  and  recognizing  that  sin  </>  =  dr! dm  gives 
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dm  dm  2r^  dm  ”  " 


(2.45) 


2.5  The  Angular  Momentum  Equation 
Multiplying  the  tangential  momentum  equation,  Eqn.  (2.11),  by  and 
rearranging  yields  a  form  of  the  angular  momentum  equation  convenient  for  use  with  the 
stagnation  enthalpy,  Eqn.  (244),  and  radial  equilibrium,  Eqn.  (2.24)  equations: 


1  S{rV,f 
2r^  dm 


(2.46) 


2.6  Closure  of  the  System  of  Equations 

In  the  preceding  sections,  seven  equations  have  been  developed  from  the  three 

component  momentum  equations  and  the  energy  equation.  These  equations  are: 

E,  viscous  body  force  Eqn.  (2.28) 

E,  viscous  body  force  Eqn.  (2.29) 

E„  viscous  body  force  Eqn.  (2.30) 

Radial  equilibrium  Eqn.  (2.24) 

Energy/entropy  Eqn.  (2.42) 

Meridional  momentum/enthalpy  Eqn.  (2.45) 

Angular  momentum  Eqn.  (2.46) 


Additional  equations  are  needed  to  close  the  system  of  equations. 

An  equation  for  density  may  be  derived  by  integrating  Gibbs’  Tds  relationship  to 


give 
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s  =  CphiT  -  R\np  (2.47) 

where  s  is  the  entropy,  R  is  the  gas  constant,  and  4  is  the  specific  heat  at  constant 
pressure.  An  arbitrary  reference  state  of  s  =  0  at  7  =  1  and  p  =  1  was  used  for 
integration  of  the  Tds  relationship.  Combining  Eqn  (2.47)  with  the  ideal  gas  equation  of 
state  yields 


(2.48) 


The  definition  of  stagnation  enthalpy  may  be  rearranged  to  give  the  temperature 
equation 


+  V? 

'  m  '0 


(2.49) 


The  streamline  slope,  0,  and  the  q-o  angle,  y,  are  determined  from  the  streamlines 
and  input  quasi-orthogonal  geometry  (see  Fig  2.4)  by 
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<|)  =  tan"*— 1 

^  "  (2.50) 

Y  =  tan"*-^] 


Through  specification  of  the  Prandtl  number,  Pr,  the  Schmidt  number,  Sc,  and  the 
mixing  coefficient,  e,  the  viscosity  and  heat  transfer  coefficient  may  be  determined  by 


H  =  Sc  p  e 

Pr 


The  mixing  coefficient  is  a  function  of  the  blade  and  annulus  geometry  and  the 
turbomachine  operating  conditions.  Those  interested  in  more  details  on  the  calculation  of  e 
are  referred  to  Gallimore  (1985)  and  Schlichting  (1979). 

The  distributed  blade  body  forces,  and  f„,  required  in  Eqns.  (2.24),  (2.45), 
and  (2.46),  depend  on  blade  geometry  for  their  calculation  (See  Appendix  A). 

For  a  given  streamline  pattern,  the  seven  equations  listed  above  combined  with 
Eqns.  (2.48)  -  (2.51)  constitute  a  complete  set  of  equations  for  streamline  curvature 
analysis.  The  input  variables  required  for  this  system  are  the  Schmidt  and  Prandtl  number, 
the  mixing  coefficient,  the  constant  pressure  specific  heat,  the  gas  constant,  and  the 
empirical  entropy  gradient  as  determined  from  blade  geometry  and  flow  conditions. 
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2.7  Inlet  Conditions 

A  consistent  set  of  initial  conditions  is  needed  to  start  the  iterative  solution  of  the 
flow  field.  The  following  assumptions  are  made  at  the  inlet  to  the  flow  annulus: 

1 .  The  inlet  station  q-o  is  assumed  to  be  radial;  therefore,  angle  y  =  0  and  dq  may 
be  replaced  by  dr  in  the  system  of  equations  (See  Fig.  2.4). 

2.  The  flow  is  inviscid  and  isentropic  prior  to  reaching  the  inlet.  As  a  result  of  this 
assumption,  Eqns.  (2.28),  (2.29),  (2.30),  (2.42),  and  (2.45)  reduce  to  zero, 
identically. 

3.  Streamlines  at  the  inlet  station  are  assumed  to  be  flat,  \.e.,  <t>  =  0,  r,  =  oo,  and 
V„  =  V,  (See  Fig.  2.4). 

4.  There  are  no  blades  at  the  inlet  station. 

5.  Stagnation  temperature  and  pressure  are  assumed  to  be  constant  across  the  inlet 
station;  therefore,  H  and  s  are  also  constant  across  the  inlet. 

These  assumptions  result  in  a  reduced  set  of  equations  for  the  inlet. 

Using  the  assumptions  listed  above,  the  radial  equilibrium  equation,  Eqn.  (2.24), 
reduces  to  the  simple  radial  equilibrium  equation 

(2.52) 

^  dr  r  dr 

Equations  (2.47)  -  (2.49)  and  (2.51),  used  to  compute  entropy,  density, 
temperature,  viscosity,  and  the  heat  transfer  coefficient,  remain  unchanged. 
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If  the  inlet  velocity  profile  contains  no  V*  component  (zero  whirl  at  the  inlet),  the 
right  hand  side  of  Eqn.  (2.52)  reduces  to  zero  and  the  axial  velocity,  must  be  constant 
across  the  inlet. 

For  the  case  of  the  inlet  velocity  profile  containing  an  arbitrary  V#  distribution 
with  radius,  Eqn.  (2.52)  must  be  integrated  across  the  inlet  to  give  the  profile  of  the  axial 
velocity,  V,. 

The  continuity  equation  is  integrated  across  the  inlet  from  streamline  to  streamline 
to  obtain  the  mass  flow  rate  which  is  used  to  satisfy  continuity  for  the  remaining  q-o%. 

The  streamline  to  streamline  integration  also  allows  the  computation  of  the  fraction  of  the 
mass  flow  which  passes  between  the  hub  and  each  streamline  outward  towards  the  outer 
casing.  The  radial  location  of  each  streamline  is  adjusted  as  necessary  to  maintain  the 
same  mass  flow  fraction  as  occurred  at  the  inlet  station.  Details  of  the  mass  flow 
integration  are  covered  in  a  following  section. 

2.8  Solution  Procedure  for  the  Governing  Equations  in  Bladed  Flow  Regions 
The  flow  annulus  of  a  turbomachii^  may  be  divided  into  unbladed  and  bladed 
regions  with  the  bladed  regions  further  divided  into  rotor  and  stator  blade  rows.  The 
solution  procedure  for  spanwise  mixing  based  on  streamline  curvature  through-flow 
analysis  for  bladed  regions  is  outlined  as  follows; 

1 .  The  velocities,  velocity  gradients,  and  streamline  pattern  from  a  previous  iteration 
or  initial  solution  estimate  are  used  to  compute  the  viscous  body  forces,  £„  E^, 
and  from  Eqns.  (2.28)  -  (2.30),  and  the  distributed  blade  body  forces,/;,^. 
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and  f„,  from  Eqns.  (A4)  and  (A5)  in  Appendix  A. 

2.  The  temperature  data  obtained  from  the  temperature  field  in  a  previous  iteration  or 
initial  solution  estimate  are  combined  with  the  empirical  entropy  gradient,  the 
velocity  gradients  and  the  streamline  pattern  to  solve  the  energy  equation,  Eqn. 
(2.42),  for  the  meridional  entropy  gradient  along  each  streamline. 

3.  The  angular  momentum  equation,  Eqn.  (2.43),  is  solved  to  give  the  meridional 
gradient  of  angular  momentum  along  each  streamline. 

4.  Using  the  gradients  determined  from  the  energy  and  angular  momentum  equations 
in  steps  2  and  3,  and  the  meridional  body  force,  determined  in  step  1,  the 
meridional  gradient  of  stagnation  enthalpy  is  computed  for  each  streamline  using 
Eqn.  (2.41). 

5.  The  entropy,  angular  momentum,  and  stagnation  enthalpy  values  are  then  updated 
for  each  streamline  and  q-o  intersection  using  the  gradients  determined  above. 

6.  The  temperatures  and  densities  are  updated  for  each  intersection  using  Eqns.  (2.48) 
and  (2.49).  The  entropy,  stagnation  enthalpy,  and  velocity  gradients  in  the  q-o 
direction  are  then  computed  using  the  data  determined  in  step  S. 

7.  The  radial  equilibrium  equation  is  integrated  to  give  the  meridional  velocity 
distribution,  V„,  along  each  q-o. 

8.  The  new  velocity  profile  computed  in  step  7  is  used  in  the  integration  of  the 
continuity  equation  at  each  q-o.  The  mass  flow  rate  thus  determined  is  used  to 
adjust  the  constant  of  integration  in  step  7  as  needed  to  achieve  the  same  annulus 
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mass  flow  as  was  determined  at  the  inlet  station.  The  solution  procedure  iterates  on 
steps  7  and  8  until  the  difference  between  the  mass  flow  rates  at  each  q-o  and  the 
inlet  lies  within  an  arbitrarily  chosen  tolerance. 

9.  The  streamline  radial  locations  are  adjusted  by  interpolating  the  mass  flow 
distributions  to  maintain  same  fraction  of  mass  flow  between  the  hub  and  each 
streamline  as  at  the  inlet  station. 

10.  Overall  convergence  of  the  solution  is  determined  by  comparing  the  streamline 
shifts  between  successive  iterations  for  the  flow  field  with  an  arbitrarily  chosen 
tolerance.  If  the  streamline  shifting  is  within  the  tolerance,  the  flow  field  solution 
has  converged.  If  not,  the  streamline  pattern  is  updated  with  the  results  of  step  9 
and  the  solution  repeated  again  starting  at  step  1 . 

2.9  Solution  Procedure  for  the  Governing  Equations  in  Unbladed  Flow  Regions 
The  solution  procedure  for  bladed  regions  differs  from  the  unbladed  regions  only 
by  the  inclusion  of  the  blade  body  forces  and  the  empirical  entropy  gradient. 
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3.  SOLUTION  METHOD,  COMPUTER  PROGRAM  DEVELOPMENT  AND 

OPERATION 

The  governing  equations  developed  in  Chapter  2  are  written  in  finite  difference 
form  and  programmed  for  iterative  solution.  The  computer  program  was  developed  for  a 
blade-free  annulus  with  inviscid  and  adiabatic  endwalls.  The  annulus  is  essentially  axial 
with  the  q-os  directed  radially  from  the  hub  to  the  outer  casing.  The  governing  equations 
are  in  the  form  derived  by  Gallimore  (1985)  which  are  a  special  case  of  the  equations 
given  in  Chapter  2. 

The  following  describes  the  main  program  involved  and  the  associated  subroutines 
which  have  been  developed.  Figure  3.1  is  a  mapping  of  the  main  program  and  the 
subroutines  called.  As  stated  previously,  the  solution  is  iterative  with  the  sequence  of 
execution  essentially  from  left  to  right  across  the  figure.  Execution  is  repeated  numerous 
times  until  convergence  is  obtained  and  Subroutine  OUTPUT  is  called  to  print  the 
solution  results. 

The  complete  FORTRAN  code  for  the  program  is  contained  in  Appendix  C. 
Comment  lines  in  the  main  program  and  subroutines  include  a  glossary  of  the  FORTRAN 
variable  names  used. 

3.1  Main  Program  and  Control  of  the  Calculations 

The  primary  purpose  of  the  main  program  is  to  control  the  order  in  which  the 
subroutines  are  called  to  carry  out  a  solution.  The  main  program  also  prints  statements 
which  list  the  subroutine  currently  being  executed  and  the  number  of  iterations  completed. 


Figure  3. 1  Main  program  and  subroutine  map 
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Subroutine  INPUT  reads  the  input  data  for  the  inlet  station  q-o  and  then  calls 
Subroutine  INLET  as  shown  in  Fig.  3.1.  Subroutine  INLET  computes  values  for  the  flow 
field  variables  such  as  static  enthalpy  and  entropy.  Subroutine  INLET  in  turn  calls 
DERIV  to  compute  required  first  and  second  derivative  values  for  flow  field  variables. 
Subroutine  VISCOUS  is  called  by  INPUT  to  compute  the  viscous  body  forces.  Eg  and  E^, 
and  the  dissipation  function,  <l>. 

Subroutine  FLOINT  calculates  the  mass  flow  rate  and  fraction  of  mass  flow  rate 
between  adjacent  streamlines  based  on  the  inlet  q-o  input  data. 

Subroutine  ESTIM  initializes  the  downstream  flow  field  variables  and  streamline 
pattern  based  on  the  inlet  q-o  input  data  and  the  hub  and  casing  geometry.  Note  that 
INPUT,  FLOINT,  and  ESTIM  are  called  only  once  at  the  start  of  program  execution. 

Subroutines  CURVE  through  RADIAL  comprise  the  iterative  portion  of  the 
program  for  the  q-os  downstream  of  the  inlet  station  and  implement  the  solution 
procedure  outlined  in  Sections  2.8  and  2.9.  A  global  iteration  consists  of  one  pass  through 
the  annulus  executing  Subroutines  VISCOUS  through  CONVERl  for  all  the  downstream 
q-o%  with  the  streamline  radial  locations  held  constant  and  is  completed  once  the  change  in 
meridional  velocity  from  one  pass  to  the  next  is  less  than  a  user  chosen  tolerance  as 
determined  by  Subroutine  CONVERl.  Subroutines  RADIAL  and  INTERP  then  compute 
the  new  streamline  locations  for  the  converged  velocity  field. 

A  universal  iteration  consists  of  executing  Subroutines  CURVE  through 
CONVER2  for  all  the  downstream  q-os.  After  a  global  iteration,  CONVER2  is  called  to 
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check  for  streamline  convergence  which  is  determined  by  comparing  the  streamline  shifts 
between  successive  universal  iterations  with  an  arbitrarily  chosen  tolerance.  If 
convergence  has  not  been  attained,  the  main  program  starts  a  new  universal  iteration  at 
Subroutine  CURVE.  Once  the  flow  field  has  converged,  Subroutine  OUTPUT  is  called  to 
print  the  converged  flow  field  values  to  an  output  data  file. 

3.2  Subroutine  Descriptions 

Each  subroutine  listing  in  Appendix  C  contains  a  glossary  of  the  FORTRAN 
variable  names  used  as  well  as  descriptions  of  the  input  required  and  the  variables 
calculated  as  output.  The  following  descriptions  contain  the  finite  difference  formulations 
used. 

Subroutine  INPUT  opens  an  existing  unformatted  data  file  called  inputl.dat  and 
reads  the  data  required  by  the  program.  Figure  3.2  shows  the  meridional  plane  with 
arbitrary  hub  and  casing  shapes  and  the  FORTRAN  variables  used  in  inputl.dat. 

Stagnation  temperature  and  pressure,  represented  by  FORTRAN  variables  TEMPIN 
and  PRESIN,  are  assumed  to  be  constant  across  the  inlet  station.  MIXCO,  the  mixing 
coefficient  described  in  Chapter  2,  is  a  function  of  blade  and  aimulus  geometry  and  the 
turbomachine  operating  conditions.  NSTRM  is  the  number  of  streamlines  in  the  flow  field 
and  NSTNS  is  the  number  of  q-as.  The  index  variable,  J,  identified  the  streamlines,  and  I 
identifies  the  quasi-oithogonals,  as  shown  in  Fig.  3.2. 

The  user  supplies  the  radial  location  of  each  streamline  at  the  inlet  station, 
RAD/US(J,1).  The  tangential  velocity  for  c*ich  streamline,  VTHETA(J,I),  and  the  axial 
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velocity,  VMID,  at  a  middle  streamline  which  is  determined  by 


JAfID  = 


+  1 


(3.1) 


The  streamline  spacing  prescribed  at  the  inlet  from  hub  to  casing  is  arbitrary. 

The  q-os  in  the  flow  field  are  limited  to  the  radial  direction,  although  the  z 
direction  spacing  is  arbitrary.  The  axial  location  of  each  q-o  is  specified  by  ZLOC(l)  and 
the  radial  location  of  the  hub  and  casing  on  the  q-o  by  RADWS(1 ,1),  and 
RADIUS(NSTRM,I).  The  arrangement  of  the  data  in  file  inputl.dat  is  explained  in 
comments  included  with  Subroutine  INPUT  in  Appendix  C. 

Subroutine  INLET  computes  the  flow  properties  for  the  inlet  q-o  using  Eqns. 
(2,47)  to  (2.49).  The  flow  is  assumed  to  be  inviscid  and  adiabatic  prior  to  reaching  the 
inlet,  and  the  streamlines  at  the  inlet  are  take  to  be  straight  and  parallel  to  the  z  axis.  For 
a  blade-free  annulus  with  radial  q-os,  the  simple  radial  equilibrium  equation  give  in  Eqn. 
(2. 52)  applies. 

Integration  of  Eqn.  (2.52)  along  the  inlet  q-o  results  in  the  following  expressions 
for  the  axial  velocity,  V,: 


+ 


(/  (3.2) 


if 
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j*l 


{'l-'j-i) 

’  1  a(rFef 

2 

dr 

J‘<  0  *■  J. 


The  integration  results  in  Eqn.  (3.3)  proceed  from  the  middle  streamline,  where  V,  = 
VM/D,  to  the  casing,  and  in  Eqn.  (3.4)  from  the  middle  streamline  to  the  hub,  stepping 
streamline  by  streamline.  The  gradient  terms  in  Eqns.  (3.3)  and  (3.4)  are  determined  in 
Subroutine  DERIV  as  explained  below. 

The  static  temperature,  entropy,  and  density  for  the  inlet  q-o  are  determined  using 
Eqns.  (2.47)  to  (2.49)  in  the  forms; 


Tj  =  T 


2c„ 


Sj  =  Cj,\n{T°)  -  R}Bip°)  ij  =  \..I^STRM) 


Py  = 


c.-R  -s, 

e 


T  R  c  R 


(3.4) 


Subroutine  DERIV  contains  the  finite  difference  schemes  needed  to  compute  the 
first  and  second  derivatives  of  the  distributions  of  flow  field  variables  in  the  axial  or 
radial  direction.  Since  the  streamline  spacing  is  irregular,  the  finite  difference  scheme 
takes  the  variable  spacing  into  account  based  on  linear  differencing  of  the  variables 
between  streamlines.  For  an  arbitrary  variable,  U(r)  along  a  q-o,  the  radial  gradient  at 
radius  r^  is  determined  with  appropriate  weighting  for  variable  streamline  spacing  as 
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dU 
dr  , 


0*1-0 

r  u  -  u. , 

J  7-1 

,  0-0-1 

\u. u.] 

7*1  7 

0*1  -  0-1 

.  0  -  0-1 , 

0*1  -  0-1 

.  0*1  -  0 . 

{j  =  2...NSTRM-l)  (3.5) 


A  one-sided  finite  difference  scheme  is  needed  to  compute  the  radial  gradients  at 
the  hub  and  outer  casing.  Applying  a  second  degree  polynomial  at  streamlines  j=l,  2, 
and  3  along  a  q-o  (Anderson,  et.  al.,  1988)  results  in  the  gradient  at  the  hub  radius  as 


dU 

u,-u,_ 

o-o' 

dTj.,  ' 

0-0 

0-0. 

.  '■3“0  . 

0-0 

(3.6) 


The  corresponding  gradient  at  the  casing  is 


dU  _  ^Nsrm-i  ~  ^Nsmt  _ 

^NSm4-\~^NSrRM 

^NSTmt-2~  ^NSntM 

^y-WSTW#  ^NSTRU-X  ~  ^NSTRIi 

^NSmi-2~^NSlK4_ 

^NSntM-2~  ^Nsmn  . 

^  ^NSTRM-l~  ^NSTltM 
^NSTIUi-2~ 


(3.7) 


The  second  derivative  of  U(r)  at  radius  rj  is  computed  as 


g‘U 

2 

I 

1 

0*1 -0-1 

.  0*1  -  0 

0  -  0-1 

a  =  2...NSTRM-l)  (3.8) 


Again,  a  one-sided  finite  difference  scheme  is  required  at  the  hub  and  casing,  with 


the  results  being 
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2 

2 

1 

» 

'•3  -'-2 

'•2-^1 

r^-r^ 

(3.9) 


and 


d^U  _  2 

^NSTRM-l  ~  ^ NSTRM 

j^NSntM  ^NSm4-2  ~  ^NSTm-\ 

^NSTRM-2  ~  ^NSTRM 

2 

^ NSTRM ^ NSTRM 

’’nSIRM-I  ~  ^NSTRM 

/nSTRM-2  ~  ^NSTRM-l, 

(3.10) 


Axial  gradients  along  a  streamline  where  needed  in  the  solution  are  computed 
using  Eqns.  (3.5)  to  (3.10),  but  with  z  replacing  r,  index  i  replacing  index y,  and  NSTNS 
replacing  NSTRM. 

Subroutine  FLOINT  computes  the  mass  flow  rate  based  on  the  input  data  at  the 
inlet  q-o.  The  integration  technique  used  is  the  trapezoidal  rule,  with  integration 
proceeding  from  streamline  to  streamline,  starting  at  the  hub  and  stepping  radially 
outward. 

The  mass  flow  rate  between  adjacent  streamlines,  j  and  j-1,  using  data  at  the  inlet 
q-o,  is  computed  by 


rh, 


P/  +  P 


{}  =2...NSTRM) 


The  total  mass  flow  rate,  ih  ,  for  the  annulus  is  the  sum  of  the  mass  flow  rates 
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lietween  the  streamlines. 

The  fraction  of  the  total  mass  flow  rate  between  the  hub  and  each  streamline  is 
then  computed  using 


^  -2...NSTRAf)  (3. 12) 

E 

*=i 


The  individual  fractions  of  the  mass  flow  rate,  therefore,  identify  the  streamlines  through 
the  annulus. 

Since  the  flow  at  the  armulus  walls  is  inviscid  and  adiabatic,  and  the  annulus  is 
blade-free,  no  global  energy  transfer  into  our  out  of  the  annulus  is  possible,  although 
energy  may  be  redistributed  within  the  annulus.  As  a  check  on  energy  conservation,  the 
mass  weighted  stagnation  enthalpy  for  the  annulus  is  computed  using 


NSTKM 


E 


m. 


m 


y-2 


(3.13) 


The  mass  weighted  stagnation  enthalpy  can  then  be  computed  for  each  downstream  q-o 
and  the  resulting  value  compared  with  the  inlet  value  as  a  check  on  the  flow  field 
calculations. 

Subroutine  ESTIM  estimates,  for  the  initial  calculations  of  the  flow  field,  the  static 
temperature,  entropy,  density,  and  streamline  radial  locations  for  the  q-os  downstream 
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trom  the  inlet  tor  i  =2  to  i  =  NSTNS,  where  NSTNS  is  the  number  of  q-os  ,  or 
computing  stations.  ESTIM  is  executed  once  only  at  the  start  of  the  calculations. 

Since  the  annulus  cross  section  may  change  from  q-o  to  q-o  as  shown  in  Fig.  3.2, 
the  streamline  radial  locations  for  the  current,  or  ith,  q-o  must  be  estimated  using 
information  from  the  inlet  q-o.  The  streamline  radial  locations  are  compute  based  on 
fractions  of  the  span  at  the  inlet  according  to 


'■y  ■  ('•.jfSTW,  - *  '•i.i  V-2  -NSTRM-1,  t-ZMSmS)  (3.14) 

^UNSTRM 


The  axial  velocity  for  the  ith  q-o  is  estimated  using  the  mass  flow  rate  determined 
at  the  inlet.  The  average  density  at  the  inlet  q-o  and  throughout  the  annulus  is  given  by 


P  = 


1 


Nsnai 


T 

NSTJiM  pi 


(3.15) 


Using  the  annulus  area  for  the  ith  q-o,  p,  and  the  total  mass  flow  rate,  m  ,  the  average 
axial  velocity  for  the  ith  q-o  is  estimated  as 


K,  =  ^ - - - —  (i  =  2...NSTNS) 

‘  P  n  {^^t,N5TKM  ~  ^  i,l) 


(3.16) 


The  stagnation  properties  are  approximated  by  simply  carrying  the  values 
downstream  from  the  inlet  station.  The  tangential  velocities  are  similarly  carried 
downstream.  The  static  temperature,  entropy,  and  density  for  the  ith  q-o  are  determined 
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using  Eqns.  (3.4) 

Subroutine  CURVE  calculates  the  streamline  curvature,  meridional  flow  angle,  and 
the  axial  and  radial  velocities,  K  and  K,  for  a  q-o.  CURVE  is  the  first  subroutine 
executed  within  a  universal  iteration  for  the  flow  field. 

The  streamline  curvature  at  the  determined  streamline  intersections  with  the  q-o  is 
computed  as 

~  (/=l..JVSmf)  (3.17) 

2 


where  the  derivatives  of  r(z)  are  determined  in  DERIV  using  values  of  the  fiinction  at 
neighboring  q-os. 

The  meridional  flow  angle  of  a  streamline,  j,  at  a  q-o  is  determined  from  the  finite 
difference  form  of  Eqn.  (2.50),  or 


♦ 


=  *  arctanf- 


-'■f-i 


\^+l 


{i^2.MSmS-l)  (3.18) 


At  the  inlet  q-o,  <^  =  0  and  for  the  last  q-o  in  the  flow  field,  4)  is  determined  using  the 
NSTNS  and  NSTNS  -  1  quasi-orthogonals. 

The  axial  and  radial  velocities  at  a  q-o  are  then  computed  as 


yj  = 

yr  = 
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(/  =  I...NSTNS) 


(3.19) 


where  V„  is  the  meridional  velocity  computed  in  subroutine  VELO. 

Subroutine  VISCOUS  computes  the  distributions  of  the  viscous  body  forces,  Eg 
and  E„  and  the  dissipation  function,  along  a  q-o.  Using  the  constraints  or  assumptions 
concerning  the  annulus  and  q-o  geometry  listed  at  the  beginning  of  this  chapter,  the 
dissipation  function  given  in  Eqn.  (2.40)  becomes 


(3.20) 


and  the  viscous  body  forces  in  Eqns.  (2.29)  and  (2.30)  are  written  as 


1  _a 

r  dr 

1  _a 

r  dr 


r\i 

r\i 


dr 


dr 


(3.21) 


The  dissipation  function  along  a  given  q-o  is  then  computed  as 


*j  ■  •*; 


fan! 

2 

+ 

--^1 

2 

(arj 

J 

\dr) 

'y, 

(j  =  IJfSTRM) 


(3.22) 
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where  the  gradient  terms  are  computed  by  Subroutine  DERIV  and  the  viscosity  is 
computed  using  Eqn.  (2.51).  The  viscosity  is  computed  in  Subroutine  ESTIM  for  the  first 
global  iteration  and  in  Subroutine  GRADIENT  for  subsequent  global  iterations. 

The  viscous  body  forces  are  computed  using 


lA 

rjdrl 


0 


dr 


_ 


O' 


j  n 


U  =  I...NSTRM) 


(3.23) 


Subroutine  GRADIENT  calls  Subroutine  DERIV  to  determine  the  meridional 
gradients  of  entropy,  angular  momentum,  and  stagnation  enthalpy  at  a  given  q-o,  and  then 
updates  these  flow  field  variables  at  that  q-o.  DERIV  is  called  again  to  compute  the  radial 
gradients  needed  for  integration  of  the  radial  equilibrium  equation  for  this  q-o  in 
subroutine  VELO. 

Under  the  constraints  or  assumptions  listed  at  the  beginning  of  this  chapter,  the 
equations  for  entropy,  angular  momentum,  and  stagnation  enthalpy,  Eqns.  (2.42),  (2.46), 
and  (2.45),  written  for  streamline  j  at  the  given  q-o  are: 
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dm  J 

•'  K  p 

J  \  ^  /j 

(j  =  LMSTRAf) 


(3.24) 


where  the  radial  gradient  terms  in  Eqns.  (3.24)  are  evaluated  using  DERIV.  The  heat 
transfer  coefficient,  k,  is  computed  according  to  Eqn.  (2.51). 

The  meridional  gradients  of  entropy,  angular  momentum,  and  stagnation  enthalpy, 
computed  by  Eqns.  (3.24),  are  used  to  update  the  flow  field  variables  at  the  current  q-o. 
The  meridional  length  of  each  streamline,  j,  between  given  adjacent  q-oSy  i-l  and  /,  is 
computed  by 

The  entropy,  angular  momentum,  and  stagnation  enthalpy  are  updated  along  the 


given  streamline  for  the  q-o  according  to 
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(3.26) 


The  change  in  angular  momentum  from  the  previous  iteration  is  relaxed  using  RF, 
an  arbitrary  relaxation  factor  chosen  by  the  user  as 

(rV.fiNEW)  =  (rV.fiOLD)  +  RF[{rV,f  -  {rV,f{OLD)]  (3-27) 

The  tangential  velocity,  V^,  is  computed  from  the  updated  angular  momentum 
above,  and,  in  turn,  the  static  temperature,  pressure,  and  density  are  updated  using  Eqn. 
(3.18). 

The  radial  gradients  of  the  updated  entropy,  angular  momentum,  and  stagnation 
enthalpy  along  the  q~o,  required  in  integration  of  the  radial  equilibrium  equation,  are  then 
computed  along  the  q-o  by  Subroutine  DERIV.  The  number  of  the  current  q-o  is  sent  to 
GRADIENT  in  the  call  list  and  the  main  program  controls  the  stepping  from  q-o  to  q-o, 
from  /  =  2  to  /  =  NSTNS. 

Subroutine  VELQ  integrates  the  radial  equilibrium  equation  to  obtain  the 
meridional  velocity  profile  along  a  given  q-o.  The  resulting  profile  is  used  to  determine 
the  mass  flow  rate  and  the  constant  of  integration  involved  is  then  adjusted  as  necessary  to 
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achieve  mass  continuity.  Once  mass  continuity  is  satisfied,  the  main  program  steps  the 
computations  to  the  next  q-o. 

Using  the  assui.iptions  of  a  blade  -  free  annulus  and  radially  directed  quasi- 
orthogonals,  the  radial  equilibrium  equation,  Eqn.  (2.24),  may  be  rearranged  to  give 


dr 


‘  B  -AVi 


(3.28) 


where 


B  = 


_  1 

dr  2r^  dr 


-  £„sin<j) 


(3.29) 


and 


^  sin<})  dV^ 

y.  dm 


(3.30) 


where  A  and  B  are  functions  of  radius,  r,  along  the  q-o. 

Equation  (3.28)  is  a  linear,  first  order,  ordinary  differential  equation  and  may  be 
integrated  by  using  a  suitable  integrating  factor. 

Multiplying  Eqn.  (3.28)  by  an  integrating  factor,  s,  and  rearranging  yields 
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2 

^sAVl  =sB 

dr 


For  the  left  -  hand  side  of  Eqn.  (3.31)  :o  be  an  exact  differential,  we  require 


dr 


+  Vm 

dr  ”  dr 


Comparing  Eqns.  (3.31)  and  (3.32),  it  can  be  seen  that 


*  =5^ 
dr 


from  which 


fAdr 

= 


s  =  e 


where  a  constant  of  integration  has  been  omitted  for  simplicity. 
Substituting  Eqn.  (3.34)  into  Eqn.  (3.31)  then  yields 


_a 

dr 


(3  31) 


(3.32) 


(3.33) 


(3.34) 


(3.35) 


Equation  (3.35)  may  now  be  integrated  along  a  q-o  from  a  point  r,  where  =  V„,  to 


obtain 
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■'*  (3  36) 

r^ 

For  a  short  distance  r,  to  r,  A  and  B  can  be  approximated  as  average  values  A  and 
B.  Therefore,  Eqn.  (3.36)  becomes 

B  *  vl^  (3-37) 

'•l 

Evaluating  the  integral  in  Eqn.  (3.37)  and  simplifying  yields 

vl  -  I  [l.O  -  (3.38) 

A 

If  S  =  0,  Eqn.  (3.38)  is  indeterminate.  Applying  L’Hopitals  rule  to  evaluate  in 
the  limit  as  S  goes  to  zero  then  gives 

(3.39) 

The  forms  of  Eqns.  (3.38)  and  (3.39)  used  to  integrate  the  radial  equilibrium 
equation  from  the  middle  streamline  to  the  tip  are 
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(J  =  J^^v-f^STRM) 


(3.40) 


and  for  integration  from  the  middle  streamline  to  the  hub 


Vl  =  ?i\i,Q  -  g-M'j-'j>i)]  +  vl  g-Mn-'j^y) 

J  A. 


(3.41) 


ij  =  Wi-1) 


where  at  is  the  meridional  velocity  (VMID)  from  the  previous  iteration  or  the 
initial  solution  estimate. 

The  change  in  the  velocity  profile  from  a  previous  local  iteration  is  relaxed  using 
the  relaxation  factor  RF  as 

vl(fiEW)  =  *  RF[K^  -  yI.(OLD)\  (3.42) 

The  integration  of  the  radial  equilibrium  equation  may  return  negative  values  for  the 
square  of  the  meridional  velocity.  If  this  occurs,  the  square  of  the  meridional  velocity  is 
set  to  1 .0 

The  mass  flow  rate  and  mass  weighted  stagnation  enthalpy  are  then  calculated 
using  Eqns.  (3.11)  to  (3.13)  as  described  for  subroutine  FLOINT.  Mass  continuity  is 
checked  by  comparing  the  difference  between  ih  (local)  and  ih  (inlet)  with  a  user  chosen 
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tolerance.  If  continuity  is  not  satisfied,  a  secant  method  routine  (Chapra  and  Canale, 
1988)  is  used  to  find  a  new  value  of  at  and  the  radial  equilibrium  equation  is  re¬ 
integrated  using  the  new  constant  of  integration.  A  local  iteration  is  designated  as  one 
integration  of  the  radial  equilibrium  equation  and  check  of  mass  continuity. 

Subroutine  RADIAL  adjusts  the  streamline  radial  locations  to  match  the  required 
fractions  of  the  mass  flow  rate  at  the  current  q-o.  The  meridional  velocities  which 
correspond  to  the  new  streamline  radii  are  then  interpolated  from  the  velocity  profile 
using  Subroutine  INTERP. 

The  change  in  radii  between  iterations  is  damped  using  a  relaxation  factor 
(Wilkinson,  1969-70).  The  optimum  relaxation  factor  as  determined  by  Wilkinson  is 


where  M.  is  the  mean  Mach  number  on  a  q-o,  span  is  the  distance  between  the  outer 
casing  and  the  hub,  and  Az  is  the  quasi-orthogonal  spacing.  Thus,  the  change  in 
streamline  radius  at  a  given  q-o  is  damped  according  to 

r{NEW)  =  r(OLD)  +  RF[t  -  r{OLD)]  (3  44) 

where  r  is  the  radius  interpolated  from  the  velocity  profile. 

Subroutine  INTERP  is  an  interpolation  routine  using  second  degree  Lagrange 
interpolating  polynomials  (Chapra  and  Canale,  1988).  XDATA  and  YDATA  are  arrays 
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of  known  values  of  NDATA  points  each,  i.e.  YDATA  =  f  (XDATA).  Subroutine 
INTERP  computes  the  ordinate  values,  YOUT,  corresponding  to  the  abscissa  points, 
XIN. 


INTERP  checks  first  to  ensure  that  all  XIN  data  points  lie  within  the  upper  and 
lower  bounds  of  XDATA  and  then  compares  each  XIN  data  point  sequentially  with  the 
XDATA  values  to  find  bracketing  values.  The  interpolation  result  is  then 


YOUT  =  fx/iv  -  XDATA{k*\ )]  \XIN  -  XDATA(k*2 )]  yPATAlk) 


[XDATA{k)-XDATA{k+\)][XDATA{k) -XDATA{k^2)] 

_ \XIN  -  XDATA{k)\  \XIN  -  XDATA{k*2)] _ 

[XDATA{k*  1 )  -  XDATA{k)]  [XDATA{k*  1 )  -  XDATA{k*2 )] 

\XIN  -  XDATAjk)]  \XIN  -  XDATAjk+l )] 
[XDATAik^l) -XDATA{k)]  [XDATA{k*2)  -  XDATA{k*  1 )] 


(3.45) 

m4rA(it+i) 

YDATA{k^2) 


Subroutines  CONVERl  and  CONVER2  check  the  convergence  of  the  meridional 
velocities  and  streamline  radii  through  the  passage.  At  each  streamline  and  q-o 
intersection,  the  change  in  velocity  and  radial  location  is  compared  with  user  chosen 
tolerances  as 


y.iOLD)  -  K 


^  tolerance  (velocity) 


(3.46) 


r(OLD)  -  r 


i  tolerance  (radial) 


If  convergence  is  not  attained,  the  old  meridional  velocity  and  radial  location  take  the 
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values  of  the  new  velocity  and  radial  location. 

The  main  program  calls  Subroutine  OUTPUT  if  both  the  meridional  velocities  and 
the  streamline  radial  locations  have  converged  or  if  the  maximum  number  of  iterations 
desired  by  the  user  has  been  exceeded.  If  the  flow  field  has  not  converged,  a  new 
iteration  is  started. 
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4.  EXAMPLE  SOLUTIONS 

Example  Problem  1  is  for  axial  flow  through  an  annular  passage.  If  no  tangential 
velocity  profile  is  included  in  the  input,  integration  of  the  radial  equilibrium  equation, 
Eqns.  (3.2)  and  (3.3),  at  the  inlet  q-o  results  in  a  constant  axial  velocity  profile.  Because 
of  this  velocity  condition  at  the  inlet,  the  viscous  body  forces,  Eqns.  (3.23)  and  the 
dissipation  function,  Eqn.  (3.20),  are  all  zero  and  the  static  temperature  profile,  Eqn. 
(3.4),  is  constant  across  the  inlet.  Thus,  the  meridional  gradients  of  entropy,  angular 
momentum,  and  stagnation  enthalpy  in  Eqns.  (3.24)  are  all  zero. 

For  constant  radius  hub  and  casing  shape,  initialization  of  the  flow  field  by 
Subroutine  ESTIM  will  result  in  carrying  the  inlet  axial  velocity  profile  and  streamline 
locations  to  the  downstream  q-os.  Integration  of  the  radial  equilibrium  equation  using 
Eqns.  (3.40)  and  (3.41)  simply  returns  the  constant  axial  velocity  since  the  streamline 
curvature  and  meridional  flow  angle  are  zero.  Because  no  changes  occur,  the  program 
will  converge  for  any  value  of  viscosity  or  convergence  tolerance  input  by  the  user. 

An  annular  passage  with  varying  hub  and  casing  shape  forces  the  estimation  of  the 
flow  field  for  q-o^  downstream  from  the  inlet  to  change.  For  a  converging  passage,  the 
downstream  axial  velocity  profiles  are  initialized  as  constant  values  from  hub  to  tip; 
however,  the  velocity  values  and  the  streamline  radial  locations  vary  from  q-o  to  q-o.  The 
viscous  body  forces,  dissipation  function,  and  meridional  gradients  of  entropy,  angular 
momentum,  and  stagnation  enthalpy  remain  zero.  Integration  of  the  radial  equilibrium 
equation  for  this  case  results  in  a  new  meridional  velocity  profile.  Thus,  the  effect  of  the 
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varying  endwall  geometry  propagates  through  the  system  of  governing  equations. 

The  program  was  run  for  the  case  of  a  converging  nozzle  with  no  tangential 
velocity  input.  The  hub  and  casing  profiles  formed  a  converging  passage  with  smoothly 
varying  first  and  second  derivatives.  The  hub  and  casing  were  symmetric  about  a 
centerline  of  the  passage.  Convergence  was  attained  after  3  global  and  2  universal 
iterations. 

The  streamline  and  quasi-orthogonal  locations  and  the  axial,  radial,  and  meridional 
velocity  profiles  from  the  converged  solution  are  shown  in  Fig.  4.1.  As  can  be  seen  in 
Fig.  4. 1  the  flow  velocities  are  not  symmetric  about  the  centerline,  although  the  passage 
is  symmetric.  This  asymmetry  occurs  because  of  the  differing  areas  between  streamlines. 
The  full  program  output  is  contained  in  Appendix  D. 

Example  Problem  2  is  for  axial  flow  with  swirl  through  a  straight  annulus.  This 
second  test  case  was  run  using  a  constant  radius  hub  and  casing  shape.  The  tangential 
velocity  profile  at  the  inlet  formed  a  free  vortex  profile.  The  program  converged  for  this 
case  after  2  global  and  2  universal  iterations.  Figure  4.2  shows  the  streamline  and  quasi- 
orthogonal  locations  and  the  axial,  radial,  meridional,  and  tangential  velocity  profiles  for 
this  test  case.  The  calculations  were  carried  out  for  21  streamlines,  however,  results  are 
shown  for  every  fifth  streamline  only.  The  program  output  is  contained  in  Appendix  D. 

Tangential  shear  stresses  on  the  radially  inward  and  outward  faces  of  the 
streamtubes  due  to  the  6  direction  velocities  should  cause  the  free  vortex  profile  to  shift 
towards  a  forced  vortex  profile.  The  tangential  velocity  field  shown  n  Fig.  4.2  shows  only 
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a  small  change  in  velocities  without  the  expected  shift  towards  a  forced  vortex  profile 
(Gallimore,  1985).  The  axial  velocity  profiles  remain  essentially  constant  throughout  the 
flow  field  with  no  streamline  radial  location  shifts. 

Cuvature  calculations  are  unstable  even  in  inviscid  flows  and  adding  viscosity 
decreases  stability  even  more.  It  is  thought  at  this  time  that  calculations  are  not  being 
carried  far  enough  downstream  for  the  flow  to  shift  to  the  forced  vortex  profile.  Work  is 
continuing  in  this  direction  at  this  time  (Grabowska,  1994). 


PROGRAM  CALCULATION  SUMMARY 
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APPENDIX  A.  DISTRIBUTED  BLADE  BODY  FORCES 

The  component  momentum  equations,  Eqns.  (2.10)  -  (2.12),  contain  the 
distributed  blade  body  force  terms ^,/j,  and^.  These  terms  are  functions  of  the  blade 
geometry  as  described  below.  Wennerstrom  (1974)  defined  the  distributed  blade  body 
force  terms  as  functions  of  both  geometry  and  the  entropy  generated  by  flow  over  the 
blades. 

Figure  A2  shows  a  point  on  an  arbitrary  blade  with  lean  in  the  intrinsic  coordinate 
system,  (n,  0,  m),  at  a  distance,  r,  from  the  axis  of  the  turbomachine.  The  quasi- 
orthogonal,  q,  lies  in  the  meridional  (m-n)  plane  at  an  angle  {4>  +  y)  from  the  n  direction 
(see  also  Fig.  2.4).  The  blade,  which  is  not  shown,  is  described  by  the  mean  camber 
surface.  The  unit  vector,  B,  is  tangent  to  the  blade  surface  at  an  angle  e  from  q,  and  lies 
in  the  (^-6)  plane.  §  is  a  unit  vector  tangent  to  the  blade  surface  at  an  angle  B  from  e„ 
and  lies  in  the  (0-m)  plane.  The  blade  body  force  due  to  pressure  differences  across  the 
pressure  and  suction  sides  of  the  blade,  fp,  is  normal  to  the  blade  surface. 

From  Fig.  A2,  the  following  unit  vectors  in  the  blade  surface  tangent  plane  may 
be  written: 


i  =  oosp  +  sinp  Sq 

&  =  cose  sin(4>  +  y)  e*,  -  sine  +  cose  oo6(4>  +  y)  e. 


(Al) 


The  unit  vector  in  the  direction  of  f,  is  found  from  the  cross  product,  B  x  §  giving 
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[ -sinp  cos€COS(<l>  +  Y)]e„  +  [cos  P  cose  cos(<j)  +  y)]^ 
[sinp  cose  sin(<t)  +  y)  +  cosP  sine] 


(A2) 


This  differs  from  Wennerstrom’s  (1974)  form  of  the  equation  due  to  an  error  in  the  cross 
product. 

The  distributed  blade  body  forces  for  the  (n,  0,  m)  coordinate  directions  are 
determined  by  separating  the  components  of  f,  to  give 


[sinp  cose  sin(<^  +  y)  +  cosP  sine] 
fe  =fp  [cosP  cose  co8(<|>  +  y)] 

4  =  ^  [sinp  cose  cos(<|)  +  y)] 


The  tangential  blade  body  force  may  be  determined  by  rearranging  the  tangential 
momentum  equation,  Eqn.  (2.11),  to  give 


(A4) 


The  magnitude  off,  can  then  be  found  from  the  tangential  component  equation  in  Eqn. 
(A3)  and  the  results  substituted  into  the  normal  and  meridional  component  equations, 
yielding 
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APPENDIX  B.  GOVERNING  EQUATIONS  IN  CYLINDRICAL  COORDINATES 

The  governing  equations  for  streamline  curvature  throughflow  analysis  (see  Section 
2.6)  may  be  written  for  cylindrical  streamsurfaces  with  streamline  coordinate,  m, 
corresponding  to  the  axial  coordinate,  z,  the  normal  coordinate,  n,  to  the  radius,  r,  and 
with  the  angle  </>  =  0  and  r^  =  oo.  Thus,  in  cylindrical  coordinates,  we  obtain. 

Viscous  Body  Forces,  (from  Eqn.  (2.13)) 
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Entropy/Energy  Equation,  (from  Eqn.  (2.41)) 
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dz 


pTF,  Sr) 
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9TV^ 


where  the  dissipation  function,  from  Eqn.  (2.40)  is  written  as 
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1  ar  r 

dr 

(B3) 


Meridional  Momentum/Enthalpy  Equation,  (from  Eqn.  (2.45)) 
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^  =  T—  +  1 

dz  dz  2r^  dz 
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Angular  Momentum,  (from  Eqn.  (2.46)) 


_L^ 
2r2  az 
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Radial  Equilibrium,  (from  Eqn,  (2.24)) 


-  T^  -  ±-  .  f 

^  dr  dr  dr  2r^  dr  ' 


(B6) 
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APPENDIX  C.  COMPUTER  PROGRAM  LISTING 
The  main  program  and  associated  subroutines  contained  in  Appendix  C  are  roughly 
in  the  order  in  which  they  are  executed  (See  Fig.  1).  Each  routine  contains  commment 
lines  which  describe  the  input  and  output  variables  required  for  operation  as  well  as 
describing  the  purpose  of  each  section  of  code. 
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RADIUS  HAS  UNITS  OP  ft  118 

RADIUS(J.l}  =  RAD1US(J, 1 ) / 12 .ODO  139  *  CALCULATE  MIDDLE  STREAMLINE  LOCATION  TO  PREPARE  FOR  INTEGRATION  OF 

140  •  THE  RADIAL  EQUILIBRIUM  EQUATION. 

WHIRL  HAS  UNITS  OP  ft^4/8ec'2  141  MID  =  INT(NSTRM/2}  *  1 

WHIRLU.l)  =  (RADIUS(J.l)  •  VTHETA  ( J .  1 )  )  •  *2  142  VM2NEWIM10I  =  VMID**2 
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74  .  the  hub  to  the  streamline,  the  last  (NSTFM)  STREAMTUBE  SHOULD  HAVE 

75  •  100%  OF  THE  FLOW  IN  IT.  MFRAC  IS  UNIT-LESS. 
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APPENDIX  D.  RESULTS  FOR  EXAMPLE  PROBLEMS 
Appendix  D  contains  the  full  program  output  for  the  two  example  problems 
described  in  Chapter  4. 
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INCREASE  IS  POSITIVE.  DECREASE  IS  NEGATIVE 
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